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Abstract 

An implicit algorithm for solving the equations of general relativistic hydrodynamics 
in conservative form in three-dimensional axi-symmetry is presented. This algorithm is a 
direct extension of the pseudo-Newtonian implicit radiative magnetohydrodynamical solver 
-IRMHD- into the general relativistic regime. 

We adopt the Boyer-Lindquist coordinates and formulate the hydrodynamical equations 
in the fixed background of a Kerr black hole. The set of equations are solved implicitly 
using the hierarchical solution scenario (HSS). The HSS is efficient, robust and enables the 
use of a variety of solution procedures that range from a purely explicit up to fully implicit 
schemes. The discretization of the HD-equations is based on the finite volume formulation 
and the defect-correction iteration strategy for recovering higher order spatial and temporal 
accuracies. Depending on the astrophysical problem, a variety of relaxation methods can 
be applied. In particular the vectorized black-white Line-Gauss-Seidel relaxation method 
is most suitable for modeling accretion flows with shocks, whereas the Approximate Fac- 
torization Method is for weakly compressible flows. 

The results of several test calculations that verify the accuracy and robustness of the algo- 
rithm are shown. Extending the algorithm to enable solving the non-ideal MHD equations 
in the general relativistic regime is the subject of our ongoing research. 
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1 Introduction 



The field of astrophysical fluid dynamics (henceforth AFD) deals with the macro- 
scopic evolution of gaseous-matter and plasmas in a wide variety of circumstances 
in astrophysics. The scope of AFD is broad, encompassing topics such as star for- 
mation, accretion phenomena onto compact and young stellar objects, dynamics 
of the interstellar medium, jets, winds and outflows emerging from around young 
stellar objects, from quasars and microquasars, supernovae explosion, y-ray bursts 
and structure formation in the universe. 

One of the ultimate aims of numerical astrophysics is to develop a black box al- 
gorithm which contains numerical solvers that are unconditionally stable, robust, 
efficient, Newtonian, relativistic and capable of treating flows that are strongly com- 
pressible, weakly incompressible, self-gravitating, radiating, magnetized multi-com- 
ponent-plasmas with high spatial and temporal accuracies on unstructured meshes 
and to provide the required solutions within the scale of hours. While this goal is 
unlikely to be achieved within the next few years, the increasing number of sophis- 
ticated numerical algorithms developed during the last two decades is remarkably 
encouraging. In particular, significant improvements have been achieved in increas- 
ing the spatial dimen sions and e nhancing the efficiency and accuracy of numerical 



algorithms (see e.g. lNagel et al.l . l2006|) . 

On the other hand, the problem of robustness of the solvers in AFD has been barely 
considered nor even seriously discussed. In this paper we discuss the robustness 
problem in AFD, present enhancement strategies and address the necessity of con- 
structing robust general relativistic implicit and radiative MHD solvers. For com- 
pleteness we review the concepts of efficiency and robustness of numerical solvers 
in computational fluid dynamics. 

A numerical solver is said to be relatively efficient if the corresponding number of 
algebraic manipulations per time step per computer-processor can be made respec- 
tively small. As a consequence, using high performance computers with large num- 
ber of computing processors, a significant progress has been achieved in improv- 
ing runtime efficiencies, provided the computing load is distributed appropriately. 
Thus, by means of modern hardware technology, the efficiency of computer-codes 
could be enhanced even without modifying the employed mathematical approach. 
This has led to uncoordinated developments of relatively large number of computer- 
codes, that may differ in efficiency and integrated physical processes but are almost 
essentially identical with respect to the mathematical solution procedure. This phe- 
nomenon can be attributed to the lack of robustness. By robustness, we mean the 
capability of a solver to be applied to a large class of problems without modifying 
the algorithmic core significantly. 

In an attempt to enhance both efficiency and robustness, the hierarchical solution 



strategy (HSS) has been suggested (Fig. Q] and [21 also see iHujeiratl (|2005|) and the 
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references therein). The HSS relies on the fact that any set of hydrodynamical, mag- 
netohydrodynamical or radiative equations are linearize-able and therefore can be 
re-written in a simple matrix equation Aq = b, where A, q, b are the coefficient ma- 
trix, vector of unknown variables and th e vector of known quantities, respect ively. 
Applying the defect-correction strategy ( Stetter . 1978 : Bohmer. Stetteii 1984 ). we 



may then re- write this equation as: A5q = d, where 6q and d denote respectively the 
vectors of small-time corrections and the defect, provided d is consistent with the 
real mathematical equations by construction. The matrix A in the latter equation 
can be replaced by a variety of matrices that correspond to a seque nce of numerica l 



approaches that ranges from purely explicit to strongly implicit (Hu jeiratl 12005). 
In this formulation, explicit methods arise as a special case, in which A is being 
replaced by the most easiest invertible-matrix: the identity matrix I. Based on this 
formulation, the Courant condition follows from the requirement that the matrix A 
should be stable-invertible. 

Therefore, strongly implicit and explicit methods are different variants of the same 
algebraic problem. While the former retain almost all the entries of the matrix A 
in the inversion procedure, the latter rely on neglecting all off-diagonal entries as 
well as crudely simplifying the diagonal elements. These methods are well-unified 
within the HSS, and that depending on the physical properties of the flow, a direc- 
tive will carefully select the entries of A that are relevant for the problem. 

In Table [T| we have summarized properties of several numerical methods. Thus, 
as long as efficiency is concerned, explicit methods are unrivalled candidates, pro- 
vided the flows are strongly time-dependent and compressible. However, due to the 
relatively large sound speed, these methods may stagnate both in modeling incom- 
pressible or even weakly compressible flows. To clarify this point, we mention that 
the time step-size in explicit methods must satisfy the Courant-Friedrichs-Lewy 
(CFL) condition: 

. . Ax Ax 
ot prn < mm } = mini 



exp U + V, v U(l + l/M) 



where the minimum function runs over all points of the domain of calculation. Here 
St exp , Ax, U, V s , At correspond to the explicit time step-size, space increment, veloc- 
ity, sound speed and the Mach number, respectively. Therefore, as At — » the flow 
becomes incompressible and the time step size approaches zero; hence a stagnation 
of the time-advancement procedure. We note that although in this case using con- 
sistent implicit solution procedures is necessary, by no means it is sufficient. Here 
it has been verified that standard implicit solvers experience serious difficulties in 
simulating low Mach number flows, typically found in the interior layers of stars, 
planets as well as around moving vehicles in the Earth atmosphere. 

The above discussion addresses the following questions: 1) Relativistic fluid mo- 
tions typically occur on the dynamical time scale. The advantages of still using im- 
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plicit solvers should be clarified. 2) Multigrid methods have been shown to display 
convergence which depends weakly on the number of unknowns in the finite space. 
In combination with nested iteration, the multigrid method can solve these prob- 
lems to truncation error accuracy in a number of operations that is proportional to 
the number of unknowns. Therefore the reason for still favouring the prolongation 
strategy over multigrid or adaptive mesh refinement needs to be explained. 3) The 
storage capacities of modern computers to date are capable of handling the entries 
of large matrices that correspond to the 3D MHD equations. Thus, the reliability 
and credibility of 3D axi-symmetric algorithms should be justified. 

In fact, there are several reasons that justify using implicit numerical procedures 
for modeling relativistic flows. In particular: 

• The set of relativistic MHD equations is generically a highly coupled non-linear 
system, which gives rise to fast growing perturbations due to non-linearities, 
thereby imposing a further restriction on the size of the time step. 

• In most of the cases the horizon of a black hole represents a geometrical singu- 
larity. The deformation of the geometry grows non-linearly when approaching 
the black hole. Thus, in order to capture flow-configurations in the vicinity of 
a black hole accurately, a non-linear distribution of the grid points is necessary, 
which, again, may destabilize explicit schemes. 

• Depending on the evolutionary conditions, non-relativistic flows may become 
ultra-relativistic or vice versa. However, almost all non-relativistic astrophysical 
flows known to date are considered to be dissipative and diffusive. Therefore, in 
order to track their time-evolution reliably, the employed numerical solver should 
be capable of treating the corresponding second order viscous terms properly. 

• The timescales of most astrophysical flows are considered to have a great dispar- 
ity. Stability requirement of conditionally stable methods however requires that 
the time step size should be a small fraction of the shortest possible timescale. 
This implies that, in order to cover a timescale of astrophysical relevance, an 
extremely large number of time steps would be required, which would give rise 
to prohibitive computational costs. Furthermore, the accumulated round off er- 
rors resulting from performing a large number of time-extrapolations for time- 
advancing a numerical hydrodynamical solution may easily cause divergence. 

• The initial conditions of most astrophysical phenomena are not known. There- 
fore, in carrying out global hydrodynamical simulations, the end solution should 
weakly depend on the initial flow configuration. Conditionally stable numerical 
methods, however, rely on time-advancing of the initial conditions. 

The latter reason may explain also why using the prolongation strategy is advan- 
tageous over classical multigrid. Worth noting is that the main building blocks of 
multigrid methods are: 

(1) Restriction, i.e., down sampling of the residual errors into coarser meshes. 

(2) Residual smoothing: reducing the high frequency errors by performing several 
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iterations, using a computationally efficient iterative procedure such as Jacobi 
or Gauss-Seidel. 

(3) Prolongation, which relies on the interpolation from the coarse onto finer 
meshes. 

The high-frequency errors here are reduced by cheap smoothing on the fine meshes, 
whereas the low-frequency errors are reduced by defect correction on the coarser 
meshes. As the bulk of the algebraic operations are made on the coarser meshes, the 
combined solution procedure is considerably efficient. However, multigrid meth- 
ods display satisfactory convergence, if the underlying flows are predominantly 
diffusion-dominated. In the case of advection-dominated flows, errors, that are re- 
sponsible for the slow convergence on the fine meshes, can be easily advected by 
the flow on the coarser meshes, thereby reducing the coarse grid correction. In the 
case of astrophysical flows, the corresponding equations may change their char- 
acter from Newtonian into ultra-relativistic or vice versa. Unlike Newtonian flows 
that are generically diffusion-dominated, relativistic flows may become predomi- 
nantly advection-dominated, depending on how large the corresponding Lorentz 
factors are. Hence, multigrid methods may fail to provide the expected rate of con- 
vergence. 

Finally we note that in order to model the formation and acceleration of relativistic 
flows in the vicinity of ultra-compact objects accurately, it is necessary to cover the 
domain of calculation by a strongly stretched mesh. Furthermore, Lorentz factors 
enhance the inner-coupling of the relativistic equations and give rise to a larger 
spectrum of non-linearities. These numerical difficulties combined with the need to 
include sophisticated magnetic and radiative processes make the construction of a 
fully 3D algorithm, at the moment, a computationally unrealizable numerical task. 

Therefore, in this paper, we do not intend to perform 3D global simulations, but 
rather focus on the algorithmic structure of unconditionally stable and robust 3D 
axi-symmetric solvers. These algorithms should enable us to search for station- 
ary or quasi- stationary solutions for the fully-coupled radiative MHD equations in 
which detailed physical processes are taken into account. 

The paper runs as follows: in Sec. 2 we describe the relativistic hydrodynamical 
equations, in Sec. 3 the transformation between the primitive and conservative 
variables is described. The numerical solution and the discretization methods are 
presented in Sections 4 and 5. In Sec. 6 we present the results of several test calcu- 
lations and end up with a summary in Sec. 7. 



2 The hydrodynamical equations in Kerr spacetime 

In the present study we intend to numerically solve the equations of hydrodynamics 
in both ultra-relativistic and purely Newtonian regimes. Unlike the usual conven- 
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Table 1 

A list of only a part of the grid-oriented codes in AFD and their range of applications. 
The matrix-equations in the first row are illustrated in Sec. 4. In these equations, q"' n+l , 
5t, A, a and d* denote the vector of variables from the old and new time levels, time step 
size, a preconditioning matrix, a switch on/off parameter and a time-modified defect vector, 
respectively, "m" in row 4 denotes the bandwidth of the corresponding matrix. 
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tion, in which the speed of light and the gravitational constant are set to unity, we 
use the sound speed as the basic measure for velocities. This is reasonable as the 
radial motion of rotating flows around a compact object can be as low as 10 5±1 the 
speed of light, whereas it is 1(T 2±1 of the sound speed. Close to the event horizon, 
all velocities become quantitatively comparable. This scaling enables the present 
algorithm to capture the structure of slow flows accurately and renders the appear- 
ance of terms that are extremely large or small due to scaling effects. Additionally, 
the present solution procedure is actually an extension of the purely Newtonian 
solver, IRMHD, into the general relativistic regime. 



2.1 The metric 



For completeness, we develop here the equations of hydrodynamics in the back- 
ground of spacetime metric of a Kerr black hole, using the Boyer-Lindquist coordi- 
nates (t, r, 9, (p). Adopting the 3+1 split of spacetime, a line element with the metric 
signature (-, +, +, +) can be written as follows: 



ds 2 = -a 2 dt 2 + h&idjj + pdtXdx* + p k dt), for j,k=l, 3. (1) 



For the Kerr metric, the line element reads: 

ds 2 = -(a 2 - p^)dt 2 + 2j3$d(pdt + h ik dx i dx k , (2) 



which corresponds to a matrix of the following entries: 



gtt g t<i 

grr 

gee 
g# 



PtP-a 2 jfy 
h„ 
h ee 
/L h M 



(3) 
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The coefficients g MV in the Boyer-Lindquist coordinates and their reiated functions 
are defined as foiiows: 



A 


= r 2 — 2r g r + a 2 


P 2 


= r 2 + a 2 sin 2 9 


r- 


= (r 2 + a 2 ) 2 - a 2 A cos 2 9 


0) 


= ?cos6> 

p 


a 


= |Va 


P 


= /3° = 0, /3? = -<f = - 


T 


= cos 2 e 

A 


h 

IL rr 


= a", h 9 e = P 2 , = 



(4) 

-2aM BH ^j 



-g = p 2 cos 9 = a VT. 



c, M B h, G, r g (= ^r 1 ), a, fi and "a" denote the speed of fight, mass of the biack 
hofe, the gravitationai constant, the gravitationai radius, the iapse and shift func- 
tions and the Kerr-spin parameter, respectiveiy. In writing these expressions, we 
made use of the coordinate transformation 9 = n/2 - 9, where we use the latitude 9 
instead of the poiar distance angie 9; hence the appearance of "cos" instead of "sin" 
in the metric terms. 



2.2 The governing equations 



Foifowing the internai energy formuf ation of lWifsonl (|1972|) and lHawfey et af.1 (|1984aUbh . 
we derive the hydrodynamicai equations from the four-veiocity u M u M = -c 2 , the 
conservation of baryonic number V^ipu 11 ) = 0, the paraifef component of the stress- 
energy conservation equation u v V fl T fiV = (to derive the internai energy equation) 
and from the transverse components (gg v + u^u^V^T^ = (to derive the momen- 
tum equations). 

For viscous flows, the stress energy tensor reads: 



= T% + {T^} =ph u»u v + Pg> AV + 



(5) 



where Tpp, T V ; S denote the stress energy tensor due to perfect and viscous flows, 
respectiveiy. P, 77, 0, are the pressure, which is calculated from the equation of 
state corresponding to polytropic or to an ideal gas, the dynamical viscosity which 
is assumed to be identical to the shear viscosity, and (= V^) which measures the 
divergence or convergence of the fluid world lines, respectively. W v = u M u v + g^ v 
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is the spatial projection tensor, whereas a corresponds to the symmetric spatial 
shear tensor: <x /JV = V^h^ + V ? mW. In the case of an ideal gas, the pressure and 
enthalpy read: 

p = (7-i)c v pr = (r-i)p6 

(6) 

h = c 2 +6 + P/p = c 2 + y&, 



where y, Cy and S denote the adiabatic index, specific heat and internal energy 
measured in the local rest frame of the fluid. For clarity, we re-write the hydrody- 
namical equations in flux conservative form: 

^ + -L|;( V=F)V r ) + -L|-( V^DV^) = (7) 
dt y=g dr y=g dQ 

dM r Id, r Id,, „ 



dt y/=gdr v ° j=gd6 



dP I M M M V dg^ 



dM e Id, r Id, „ 

dt y=g dr y=g d# 



d/> 1 ( M^ } ^ + L2 , (9) 

1 Jiff j +jLZ "" W 



d<9 2^ M< d6 

dM A 1 d , , _ 1 d 



+ V^M^) + -=— ( V^M^) = L2£ s (10) 



dt ^dr v ° 9 ^d6 



+ + T-A, (11) 



where D = pw' is the modified relativistic mass density. are the four-momenta: 
(M t , M r , Me, M^) = D(u t , u r , u#, u^), where D = Dh, and u' is the time-like velocity, 
V M = u^/u' is the transport velocity. L2f,- S are the spatial projections of the viscous 
stress energy tensor T V j s (see Eq. [5]) in the respective direction. These are obtained 
from the projection of the viscous tensor along the vector normal to the hyperspace, 
i.e., constant in time: 

^vis — Vfi 1 Vis _ U M 1 Vis + 1 fiA 1 Vis' 
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where £ = {r, 9, ip\. corresponds to the spatial divergence of a tensor taken in the 
Boyer-Lindquist coordinates and Y*^ are the Christoffel's symbols of the second 
kind. 

From the collection of the numerous viscous terms, we only consider the dominant 
second order operators, that are set to degenerate into the classical non-relativistic 
formulation of the Navier- Stokes e quations if the sound speed becomes smaller 
than a critical value (|Tassoul . 1978 ). 



The viscosity coefficient here is based on the alpha-turbulent description, a t , mod- 
ified to respect causality. Hence the dynamical viscosity reads: 



D D 

T] = p<l>< Vturb >=p<l> [aVJu'] = —^[a <1>V S ] = —^v t , (12) 

(u ) (U ) 



where <> denotes a turbulent mean, v tLir b is the relativistic turbulent velocity and v, 
the relativistic turbulent velocity coefficient and a, is a constant of order unity. 



Equation (fTTT) describes the time-evolution of the relativistically modified internal 
energy fi d = D C v T, where T is the plasma temperature. <J>, Y, A correspond 
to heat function due to turbulent dissipation, other heating and cooling functions, 
respectively. Using the tran sformation q /; = R IIV Q V , we may define the transport 
velocities Y M as follows (see lHawley et all Il984allbh : 



u r = 

Ug = 
U t 

II ■ = 



8rr W 

gee u 6 

g tt u< + g ttp 

gt<p u l + u 4 



= g rr u< {Vic) 

= gee u< (V e /c) 

= u t [g tt + g tl/> (Vl>/c)] 

= u'ig^+g^iVf/c)]. 



(13) 



The corresponding relativistic 4-momenta then read: 



M' 


= Du* = D 


M r 


= D[V7c] 


M° 


= D[V e /c] 




= D[V*/c], 



(14) 
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from which the covariant 4-momenta may be obtained: 



M t 



g„ M t + g t<p M* 



M r 



grr M r 

geeM 6 

g H Mt + g4>t M t . 



(15) 



< 



M, 



We note that by using this formulation of the HD-equations in combination with fi- 
nite volume discretization, mass and momenta are conserved up to small discretiza- 
tion errors. This is necessary in order to assure that inflowing non-rotating matter 
would not gain angular momenta though it must rotate in the ergosphere. 



2.3 Non-dimensional formulation 

The algorithm presented here should be capable of modeling the time-evolution 
of hydrodynamical flows both in the non-relativistic as well as in the extreme- 
relativistic regimes. In order to avoid the appearing of extremely small coefficients 
in the equations, the scaling variables listed in Table © are adopted. 



Scaling variables 




Example (supermassive BH) 


Mass: 


M 


3 x 10 8 M Q 


Accretion rate: 


M 




Distance: 


R 


3R S 


Temperature: 


t 


10 8 K 


Velocities: 


V 


V S 


Density: 


P 


1(T 12 gcnT 3 



Scaling variables that might be used to re-write the equations in non-dimensional form. In 
this table R$ (= 2r g ) stands for the Schwarzschild radius (r g is the gravitational radius) and 
Vs is the sound speed. 

We now introduce the following additional non-dimensional parameters: 



Table 2 




(16) 



where a is the black hole spin. 
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The normalization of the 4-velocity and momentum yields: 

g„ + 2 (|) 2 ^ V* + 6 6 V r( n 2 + ^(V ) 2 + (|) 2 ^(V*) 2 



(t/ ( ) 2 



= -1,(17) 



and 



g « M f + 2 (|) g^MtMf + ee 2 g rr Ml + £f>u\ + (|) «*Mj = -D 2 , (18) 



where 



V 


= -l/a 2 




- ar/(aS) 2 




= [1 -(£)(§?)] 




- l/g rr = A/p 2 




- I/gee - 1/P 2 . 



(19) 



We may write the equations of hydrodynamics in non-dimensional formulation in a 
manner that they smoothly adopt the Newtonian form in the non-relativistic regime: 
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-<? fi d V) 
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u l V 



(20) 



(21) 



(22) 
(23) 
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+ o + r - a, 



(24) 



where U = V r ^ V = Y e y/gTe, V = y^, h = 1 + (^)T. 
For flows approaching rotating black holes, the angular momentum is defined as: 
= Dyg^(V^ + Vp DE ), where V FDE denotes the rotation of the flow that is 
induced due to the frame-dragging effect: V EDE = J3^/ -yjg^- 

Note that the radial velocity in this formulation approaches the speed of light as the 
matter crosses the event horizon. 



3 The primitive variables 



The above set of equations describes the time-evolution of the conserved quantities 
D, M t and S d . However, the equation of state, the rate of transport, the applied work, 
cooling and heating rates are functions of essentially the primitive variables p, V 
and T. 

Since the relation between the primitive and conservative variables is rather non- 
linear, an iterative solution procedure should be employed. 

We note that the 4-momenta must satisfy the normalization condition: M^M^ = 
—M — -D . This is equivalent to solve the following equation for M t : 

MpM? = g^M^My = g n M 2 + 2g#M t Mf + g rr M 2 r + g ee M 2 e + g^Mj = -D 2 . (25) 



Taking into account that the quantities M r , M 6 , M<p are known at the end of each 
time step, we may substitute them in Eq. ([24]) to obtain a quadratic equation for M t , 
i.e., 

KM? + 3M t + 3 = 0, (26) 

where N = g", 2 = 2g^M^ and 3 = g rr M 2 r + g ee M 2 + g^M 2 . 
Having obtained M t , the contravariant quantity M' can be computed using the trans- 
formation: M' = g"M t + g'^Mp, whereas the global Lorentz factor is obtained from: 
u' = M'/Dh. Using Equation (1241) . the pressure can be calculated then from the 
relation: P = (y - l)6 d /u\ 



4 The hierarchical solution strategy - HSS 



The set of hydrodynamical equati ons are solved within the context of the hierar 



chical solution strategy (HSS, see iHujeiratl (|2005|) ). HSS is based on constructing 
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Fig. 1. A schematic description of the hierarchical clustering of the coefficient matrix em- 
ployed by the algorithm. Here, a cluster of coefficients is computed in the very early stage. 
The matrix-generator then selects the entries to be used for constructing the matrix co- 
efficient appropriate for the solution procedure. Depending on the matrix used, the solu- 
tion method may range from purely explicit to fully implicit. Interchange between solution 
methods is possible, as modifying, adding or removing entries is directly maintainable. 

a coefficient matrix A, which results from linearizing the complete set of equations 
in a fully implicit manner. Noting that the conservative formulation of the HD- 
equations yields a matrix coefficient that is highly sparse, it is reasonable to design 
a procedure which selects the entries for constructing the approximate matrix A 
most appropriate for the flow problem. Depending on the structure of A, a suitable 
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Fig. 2. A schematic description of the hierarchical structure of the algorithm. Stage I cor- 
responds to the implicit operator splitting approach (IOS), which is most appropriate for 
following the early time-dependent phases of the flow. The solutions obtained can then be 
used as initial conditions for Stage II, where the hydro-equations are solved as a single 
coupled system, followed by the magneto component, which is again solved as a single 
coupled system. Here, high spatial and temporal accuracies in combination with the pro- 
longation/restriction strategy may be used. Similarly, the solutions obtained in this stage 
may be used as starting solutions for Stage III, where steady solutions for the fully coupled 
set of equations consisting of the zero moment of the radiation field and the MHD equations 
are sought. In this stage, pre-conditioned Krylov sub-space iterative methods are consid- 
ered to be robust and efficient. The very last stage, Stage IV, corresponds to the case where 
solutions for the internal energy equations weakly coupled with the 5D radiative transfer 
equation are sought. 
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iterative method within the context of defect-correction method may be employed 
to assure consistency and convergence. 



For example, if we want to simulate a two-dimensional weakly compressible, non- 
magnetized and non-radiating flow between two concentric spheres, then the above- 
mentioned procedure is set to select the entries from the cluster of coefficients that 
corresponding just to the equations to be solved (see Fig. [B, which are used then 
to construct the preconditioning A. 

To clarify the procedure, we re-write the set of equations in a conservative vector 
form: 

dq -* -> -> 

-1 + L r , n .F + L Bfie G = f, (27) 



where F and G are fluxes of q, and L r iT , Le^e are first and second order operators 
that describe the advection and diffusion of the vector variables q in r and 6 direc- 
tions. / corresponds to the vector of source functions. 



In order to enhance the mathematical consistency and increase the spatial and tem- 
poral accuracies of the numerical scheme without a substantial increase of the com- 
putational costs, we adopt the so-called prediction-correction iteration procedure. 
Therefore, we re-write Eq. (12"7T ) in residual form: R = ^ + L VU F + L eM G - f = 0, 
and adopting a five star staggered grid discretization, we may apply the Newton- 
linearization to calculate the Jacobian, / m i,ni - where ml, nl are integers that 
run over the number of equations and variables. The solution can be obtained then 
as follows: 

<7 - 4 - J ml,nl K ' 



where i is the iteration level. By inspection of the Jacobian J, it can be easily verified 
that it has the following block matrix structure: 

+ «Vi + D< Hi + = RHS °f < 28 > 



where Sq = q l+l - q', and which, in the linear case, reduces to time-difference of 
q. The subscripts "j" and "k" denote the grid-numbering in the r and 6 directions, 
respectively, and RHS n = [f - L Tn F - L 0go G] n . SJ ,e and S r '° mark the sub-diagonal 
and super-diagonal block matrices, respectively. D T ' e corresponds to diagonal block 
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matrices. 

To outline the directional dependence of the block matrices, we re- write Eq. (1271) 
in a more compact form: 

S <%k+i 

+S_ T 6q hhk +D mod 6q iM +Tdq^ = RHSf k (29) 
+S_°6q hk _ 1 , 



where D raod = I /St + D 1 + D e . 

This equation gives rise to at least four different types of solution procedures: 

(1) Classical explicit methods are very special cases in which the sub- and super- 
diagonal block matrices together with D T and D e are neglected. The only ma- 
trix to be retained here is (l/6t) x (the identity matrix), i.e., the first term on 
the LHS of Eq. (1281) . This yields the vector equation (see M5/Fig.[T]): 

[|]<% k = RHSf k . (30) 

(2) Semi-explicit methods are obtained by preserving the diagonal entries, dj k , 
of the block diagonal matrix D mod (see M4/Fig. [T]). This method has been 
verified to be numerically stable even when large Courant-Friedrichs-Levy 
(CFL) numbers are used. In particular, this method is absolutely stable if the 
flow is viscous-dominated. 

(3) Semi-implicit methods are recovered when neglecting the sub- and super- 
diagonal block matrices only, but retaining the block diagonal matrices (see 
M3/Fig.[T])- In this case the matrix equation reads: 

D mod 6q jM =RHSl k . (31) 

We note that inverting D mod is a straightforward procedure, which can be 
maintained analytically or numerically. 

(4) A fully implicit solution procedure requires retaining all the block matrices on 
the LHS of Eq. (1281) . This yields a global matrix that is highly sparse (Ml/Fig. 
[[])• In this case, s emi-direct methods suc h as the "Approximate Factorization 



Method" (AFM: Beam . Warming, 1978) and the "Line Gauss-Seidel Relax 



ation Method" (-LGS: Ma cCormackl . Il989|) are considered to be efficient 



preconditioners for solving the set of radiative MHD-equation s within the 
context of defect-correction iteration method (see lHujeirat ( 2005 ) and the ref- 



erences therein). Furthermore, Krylov sub-iterative methods may prove to be 
more efficient and robust than the above-mentioned semi-direct methods. 

In the case that only stationary solutions are sought, convergence to steady state can 



be accelerated by adopting the so-called "Residual Smoothing Method" (Hujeirat, 
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2005) This method is based on associating a time step size with the local CFL- 



number at each grid point. While this strategy is efficient at providing quasi- stationary 
solutions within a reasonable number of iterations, it is incapable at providing phys- 
ically meaningful time scales for features that possess quasi-stationary behaviour. 
Here we suggest to use the obtained quasi- stationary solutions as initial configu- 
ration and re-start the calculations using a uniform and physically relevant time 
step. 



5 Numerical techniques 



In this section we briefly describe several algorithmic steps for solving the continu- 
ity equation and the generalization procedure. 



V, 



],k+l 



J* 



--a 



k+l 



Uj+l,k 



J x j+1 

Fig. 3. The location of the variables using the staggered grid discretization. 



(1) The continuity equation is discretized using the staggered grid strategy within 
the context of finite volume philosophy (Fig. [3]). 



y/-g/ 8rr DU 



1 d 



[J yFgWj 



O+i 



,(32) 



1 d 



( \/-glgee) 



DV 



J0* 



(33) 
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where 



D J,k+fj,k if ' M < « , 



D j* + fj,k if V'm<0 
i»M i if > 



(34) 



The functions f and / e are corrections for maintaining higher order spatial 
accuracies. 

(2) In order to obtain second order temporal accuracy, we write the continuity 
equation as follows: 



6D 
~~8t 



+ &Ll n+l (D) + (1 - &)Ll"(D) = 0, 



(35) 



where & denotes the parameter of the stabilized Crank-Nicolson method for 
achieving second order temporal accuracy. Ll h resembles the advection oper- 
ator at the new time level (n+1) and the old time level (n) and 6D = D n+1 -D n . 
Taylor-expanding the variable D n+i in time and considering first order terms 
only, the continuity equation gets the following form: 



en 

— + &Ll n+1 (6D) + L\ n (D) = 0, 
St 

(3) Define the defect dj ;k at every grid point: 



6D 



(36) 



(37) 



where the subscript "high" means that the transport operators are evaluated 
using a spatially accurate advection scheme. 



(4) Define at each grid point the following operator: 

LD = ~- VLV^(D) - (1 - ^)Ll^ gh (D)|^. 

Compute the following entries at each grid point: 

_ OLD _ OLD r _ OLD 

a it - , JJj,k - ;) r , , * & - 



(38) 



- j ' k azw 



OLD 



(39) 



jfc+1 



(5) In the one-dimensional case, the following matrix equation should be solved 
at each grid point: 



\SD j+hk ) 



, where j=l, J, and k= const. (40) 



\dj+i,k J 
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For J number of grid points in the radial direction, this yields the tri-diagonal 
matrix equation: 



S_ 2 £> 2 S 2 



V 



6Di 

6D 2 



d 2 



for constant k 



(41) 



S_j £>j J \sd n J yd N J 



Although this matrix equation is linear in D, it should be solved iteratively to 
recover the high spatial accuracy on the right hand side. 

Similarly, if the continuity and the radial momentum equation are to be 
solved in one dimension as a coupled system, we may obtain the following 
relation at each grid point: 





LI 1,-1 


L12 M 


LUj 


LUj 




LU j+l 


) 




£21,-1 


L22 h x 


L2\j 


L22 s 


L21 ;+1 


L22j+\ 






( 



d M 

d M 
d u 
d M 



\ 



j-1 



(42) 



for j = 1, J and k = const., 

OLm I 



where Lmni = ^1,=/. Specifically, LI is the density equation and L2 the 
momentum equation. 

For J number of points this yields a tri-diagonal block matrix, in which each 
block has the dimension 2x2. 

For a given set of equations in one dimension, we have just to replace the 
above 2x2 block matrix by a square block matrix whose dimensions are 
Neq x iVeq, where A^ eq is the number of unknown variables: 



/ 



d 2 



(43) 



q here is a vector of A^ eq entries. 
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The extension into two-dimensions gives rise to a matrix equation of the 
following form: 



dj.i 

di,2 



(44) 



\^J,K J \dj,K- J 

The matrix A has a similar structure as Ml in Fig. \T\ This matrix equation is 
solved iteratively, using a non-direct inversion procedure. 



6 Test calculations 



The verification tests of the Newtoni an version of the present algorithm have been 



presented in a series of papers (see H ujeiratl ((2005) and the references therein). 



Nevertheless the modifications made here are serious and deserve appropriate test 
calculations to ensure bug-free runs as well as a consistent re-production of the 
results in the Newtonian regime. 

In the following we briefly mention several of the test calculations performed: 

• The shock tube problem - STP 

In the case of low fluid- velocities, the modification made should enable capturing 
of shocks propagating at sub-relativistic speeds, irrespective of the accuracies 
use d. Therefore, we have applied the algorithm to the well-known Sod-problem 
(see Hujeirat ( 1995|) and the references therein). Fig. |4] shows that the algorithm 



is indeed capable of re-producing Sod's solution with high accuracy. 
The ultra-relativistic shock tube problem 

The speed of the shocks in the Sod's problem can be made arbitrary large, de- 
pending on the initial ratio of the pressure in the tube. While non-relativistic 
solvers may produce propagating velocities that exceeds the speed of light, a 
conservative and accurate relativistic solver should produce velocities that can 
be extremely close to but never exceed the speed of light. 
In Fig. [5] the one-dimensional profiles of the density, velocity, temperature and 
Lorentz factor u l are displayed. These profiles ag ree qualitatively with the ana- 



lytical solution of the relativistic STP provided by Ma rti and Mullen (|2003i) . In a 
forthcoming paper, we intend to quantitatively compare the profiles for extremely 
large Lorentz factors. 

Fig. \5\ demonstrates the strong robustness of the algorithm and its capabil- 
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ity to capture the propagation of extreme ultra-relativistic shocks in which the 
Lorentz factor is of order 1000. Such robustness is essential to enable modeling 
jetted Gamma-Ray bursts, where the Lorentz factors are in the excess of several 
hundreds. 



Non — Relativistic Shock Tube Probler 
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Fig. 4. The classical non-relativistic shock tube problem. The profiles of the density, ve- 
locity, temperature and pressure are displayed. The advection scheme used here is second 
order in time and third order in space. 200 uniformly distributed finite volume cells are 
used. 

• Relativistic Bondi accretion onto Schwarzschild black holes 

This problem is appropriate for testing the capability of the solver at treating 
transonic stationary accretion flows onto Schwarzschild black holes, assuming 
perfect spherical symmetry. This problem has been investigated by several au- 
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Fig. 5. The ultra-relativistic shock tube problem. The radial distributions of the density, 
radial velocity, temperature and the modified Lorentz factor U' are shown. The accuracy of 
the scheme and the number of points are identical to those in Fig. |4] This calculation shows 
that shocks propagating with Lorentz factors of order 1000 can be safely treated with our 
algorithm. 



thors (see iMichel ( 1972 ). also see Hawley et al. ( 1984a b) for a comprehensive 
description of the numerical treatment). In this problem, a constant flux of an 
ideal gas is said to be accreted by a non-rotating black hole. Depending on lo- 
cation of the outer boundary and on the temperature of the flow, the initially 
subsonic inflow should make a transition into the super-sonic regime at a spe- 
cific radius, which appears to be determined entirely by the constant of motion. 
On the other hand, the Lorentz factor of the flow as it crosses the inner boundary 
should approach the speed of light, depending on how close the inner radius is 
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to the event horizon. In Fig. [6] we display the radial distributions of the velocity, 
density, temperature, Lorentz factor and the Mach number, which clearly well- 
agree with the known analytical solutions. In obtaining these results we used a 
pseudo time-stepping scheme to enhance convergence. The very last time step 
size in this calculations corresponds to Courant number 2000, approximately. 
• Standing shocks around black holes 

The purpose of this test is mainly to examine the capability of the algorithm at re- 
producing the formation of the two-dimensional curved standing shocks around 
a Schwarzschild black hole that have been obtained using the Newtonian version 
of the algorithm. 

This problem is similar to the forward facing step in computational fluid dy- 
namics. Here a cold and dense disk has been placed in the innermost equatorial 
region: [1 < r < 10] X [-0.3 < 6 < 0.3] (see Figures |7] and [8]). Vanishing in- 
and out-flow conditions have been imposed at the boundaries of the cold disk. 
The gas surrounding the disk is taken to be inviscid, thin, hot and non-rotating. 
The cold disk here serves as a two-dimensional barrier that disturbs the gas from 
otherwise a spherically symmetric freely falling flow onto a Schwarzschild black 
hole and, instead, it forces the inflow to form a curved shock which eclipses the 
cold disk. 

In solving the HD-equations, an advection scheme of third order spatial accu- 
racy and of first order accuracy in time has been used. 

Hence the scheme is taken to be highly diffusive in time in order to damp os- 
cillations and to accelerate convergence into steady-state. The domain of calcu- 
lation is sub-divided into 200 strongly-stretched finite volume cells in the radial 
and 60 in the horizontal direction. In Fig. [F]the ID radial and horizontal pro- 
files, the 2D configuration of the density, temperature and the velocity field are 
shown. Indeed, the algorithm shows that it is numerically stable and capable of 
capturing steady-state shocks with complicated shock structures even for large 
CFL-numbers. 



7 Summary 

In this paper we have extended the previous Newtonian implicit algorithm to enable 
solving the hydrodynamical equations in general relativity. The 3D axi-symmetric 
hydrodynamical equations have been presented in the background of a Kerr metric 
of a black hole using the Boyer-Lindquist coordinates. The equations have been 
formulated in conservative form and subsequently solved numerically, using the 
finite volume formulation. The new extension can be well accommodated within 
the hierarchical solution scenario, in which the degree of implicitness can be made 
dynamical, depending on the hydrodynamical problem in hand. In particular, for 
modeling strongly time-dependent astrophysical flows, such as moving shocks, the 
pre-conditioners used are tri-diagonal matrices that are solved successively. Al- 
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though the computational costs per time step may be one order of magnitude larger 
than their explicit counterparts, this can be compensated through a reduction of the 
overall number of time steps required to recover a physically reliable time scale. 

On the other hand, the efficiency and robustness of the HSS are superior, if the 
solutions sought are stationary or quasi-stationary, irrespective of whether the flow 
is dissipative or not. 

Finally, a unification scheme for various numerical methods has been presented. In 
particular, the HSS algorithm enables the construction of a large variety of solvers, 
in which the degree of implicitness may range from purely explicit up to strongly 
implicit, depending on the physical properties of the underlying flow problem. 
Thus, the HSS is actually a unified algorithm for treating weakly compressible, in- 
compressible, time-dependent, time-independent, radiative, magnetized non-dissi- 
pative or strongly dissipative flows. As a consequence, using the HSS algorithm, 
we are able to save a large number of working hours which otherwise would go in 
designing different solvers for different physical problems. 

In a subsequent paper, we intend to discuss and describe the inclusion of the mag- 
netohydrodynamical equations in general relativity into the present solver. 
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Fig. 6. The Bondi accretion problem onto a Schwarzschild black hole. The profiles of the 
radial velocity, the relativistic density, temperature, the Lorentz factor U' and the radial 
Mach number. In the lower-right panel profiles of the Courant number (solid line) and the 
corresponding residual (dashed line) versus the iteration number are displayed. Although 
the problem is spherically symmetric, the calculations have been carried out using 200 
grid points in the radial and 30 in the horizontal-direction. The accuracy of the advection 
scheme is set to be first order in time and third order in space. This test enables us to 
examine the capability of the algorithm at capturing steady solutions that are essentially 
one-dimensional using a 2D numerical scheme. We have verified that the 30-profiles in the 
radial direction obtained at different are identical to machine- accuracy. 



28 




-10 12 3 4 5 

Radius (In R,„ unilal 



Fig. 7. The 2D distribution of the temperature (in units of 10 K) of a freely falling non-rela- 
tivistic gas onto a Schwarzschild black hole surrounded by a static cold disk (top panel). In 
this figure, color gradients run as follows: red color corresponds to large temperature- val- 
ues, green to intermediate and blue to low values. The distribution in the second and third 
panels have been obtained using the general relativistic version of the algorithm. Here the 
inflowing matter across the outer boundary has the temperatures 10 9 K (middle) and 10 10 K 
(bottom). 
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Fig. 8. Distribution of the grid points used in the calculations. A strong non-uniform dis- 
tribution of the grid points has been constructed to enable accurate capturing of standing 
shocks surrounding the cold disk. The tensor-product mesh consists of 275 finite volume 
cells in the radial and 130 in the horizontal direction, respectively. In the lower panel the 
profiles of the density, temperature, radial and vertical velocities along the equator and 
horizontal along the constant radius r = 2.5 are shown. 
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